NetNews Offline 2/NetNews Offline Volume 2.iso/news/comp/lang/c-part1/3395

Internet Message Format | 1996-08-05 | 6.3 KB

From: [e-mail] (MARNEFFE Thierry G.)
Subject: RE: PI Algorithm...
Date: 28 Jan 96 14:14:12 -0800
References: <4ebmpu$[e-mail]>
Message-ID: <[e-mail]>
Path: news.msn.com!msn.com
Newsgroups: comp.lang.c
Organization: The Microsoft Network (msn.com)


PI Computing: here is a solution ( up to 64.000 digits)

[e-mail]
------------------------

1. Mathematical Background
--------------------------

   Pi can be computed easily with the Stormer Formula:

     Pi = 24 arctg 1/8 + 8 arctg 1/57 + 4 arctg 1/239

   or the Gauss Formula:

     Pi = 48 arctg 1/18 + 32 arctg 1/57 - 20 arctg 1/239

   with
                  1     1     1     1     1
     arctg 1/x = --- - --- + --- - --- + --- ...
                  X    3X3   5X5   7X7   9X9

   Considering the deinition of Arctg function, the Gauss formula
   is slightly faster than the Stormer Formula and is choosen for
   computation.

   Here is what you can expect with a 486 100 Mhz:

           1000 digits     1 seconds
           2000 digits     4 seconds
           4000 digits    14 seconds


2. Definitions and Declarations
-------------------------------

  A. Data

    // Number of Decimals that must be computed ( multiple of 5)
    long      m_lDecNbr;

    // Number of Terms to compute in the Gauss Formula to obtain the
    //        required number of exact digits
    long      m_lTermNbr;

    // Pointer to Pi Digit Array
    long*     m_plPi;

  B. Functions

    // Allocate Memory for m_plPi Array
    //   Parameter: Ouput of GetArrayLen Function
    BOOL  AllocPiMem( long lLen);

    // Free m_plPi Array
    void  FreePiArray();

    // Compute Pi (...)
    void  ComputePi();

    // Get Length of Arrays for Storage of Pi and Temporary Digits
    long  GetArrayLen()
           { return (( m_lDecNbr / 4L)  + 5L); }


3. Source Code
--------------

void ComputePi()
{

    // Utility Members

    long  lLen, lEvalTerm, lTerm, lDivider, lRemainder;
    long  lComputedTerm, lEndTerm, lTermSum, lFinalSum, lFinalTerm;
    long  lStartTerm1, lStartTerm2, lStartTerm3, lDiv1, lDiv2, lDiv3;

    // Flag for Operation Type ( Addition or Soustraction)

    BOOL  bAdd;

    // Array for Storage of Intermediate Computations
    long* plArctg1;
    long* plArctg2;
    long* plArctg3;


    // Alloc Array for Pi Storage

    lLen = GetArrayLen();

    if ( !AllocPiMem( lLen))
        return;

    // Allocate Memory for Intermediary Data

    plArctg1 = new long[lLen];
    plArctg2 = new long[lLen];
    plArctg3 = new long[lLen];

    // Init Array

    for ( long i = 0L; i < lLen; i++)
    {
        plArctg1[i] = 0L;
        plArctg2[i] = 0L;
        plArctg3[i] = 0L;
    }

    bAdd = TRUE;

    //Set Starting Value for Terms Computing

    plArctg1[0] = 864L;        /* 48  * 18  */
    plArctg2[0] = 1824L;       /* 32  * 57  */
    plArctg3[0] = 4780L;       /* 20  * 239 */

    // Set Dividers

    lDiv1 = 324L;              /*  18 * 18  */
    lDiv2 = 3249L;             /*  57 * 57  */
    lDiv3 = 57121L;            /* 239 * 239 */

    // Set Start and End Term

    lStartTerm1 = 0L;
    lStartTerm2 = 0L;
    lStartTerm3 = 0L;

    lEndTerm = m_lTermNbr;

    lEvalTerm = ( m_lDecNbr / 4L) + 5L;

    // Start Computing

    lDivider = 1L;

    lComputedTerm = 0L;

    while ( lComputedTerm < lEndTerm)
    {

        // Compute Series for Term Arctg 1/18

        lRemainder = 0L;

        for ( lTerm = lStartTerm1; lTerm < lEvalTerm; lTerm++)
        {
            lTermSum = ( lRemainder * 10000L) + plArctg1[lTerm];

            plArctg1[lTerm] = lTermSum / lDiv1;

            lRemainder = lTermSum - ( plArctg1[lTerm] * lDiv1);
        }

        // Compute Series for Term Arctg 1/57

        lRemainder = 0L;

        for ( lTerm = lStartTerm2; lTerm < lEvalTerm; lTerm++)
        {
            lTermSum = ( lRemainder * 10000L) + plArctg2[lTerm];

            plArctg2[lTerm] = lTermSum / lDiv2;

            lRemainder = lTermSum - ( plArctg2[lTerm] * lDiv2);
        }

        if ( lStartTerm2 < lEvalTerm)
            if ( plArctg2[lStartTerm2] == 0L)
                lStartTerm2 += 1L;

        // Compute Series for Term Arctg 1/239

        lRemainder = 0L;

        for ( lTerm = lStartTerm3; lTerm < lEvalTerm; lTerm++)
        {
            lTermSum = ( lRemainder * 10000L) + plArctg3[lTerm];

            plArctg3[lTerm] = lTermSum / lDiv3;

            lRemainder = lTermSum - ( plArctg3[lTerm] * lDiv3);
        }

        if ( lStartTerm3 < lEvalTerm)
            if ( plArctg3[lStartTerm3] == 0L)
                lStartTerm3 += 1L;

        // Compute Sum of Terms

        lRemainder = 0L;

        for ( lTerm = lStartTerm1; lTerm < lEvalTerm; lTerm++)
        {
            // Add Arctg series

            if ( bAdd)
               lTermSum = plArctg1[lTerm] + plArctg2[lTerm] - plArctg3[lTerm];
            else
               lTermSum = plArctg3[lTerm] - plArctg1[lTerm] - plArctg2[lTerm];

            lFinalSum = ( lRemainder * 10000L) + lTermSum;

            lFinalTerm = lFinalSum / lDivider;

            lRemainder = lFinalSum - ( lFinalTerm * lDivider);

            // Add Series Sum to Current Value

            m_plPi[lTerm] = m_plPi[lTerm] + lFinalTerm;

            // Check for Carry Over

            if ( m_plPi[lTerm] < 0L)
            {
                m_plPi[lTerm] += 10000L;
                m_plPi[lTerm - 1L] -= 1L;
            }
            else if ( m_plPi[lTerm] > 10000L)
            {
                m_plPi[lTerm] -= 10000L;
                m_plPi[lTerm - 1L] += 1L;
            }
        }

        // Update Index of First Term to compute if required
        // This allow substantial reduction of computing time as this
        // avoid working on terms equal to 0L

        if ( lStartTerm1 < lEvalTerm)
            if ( plArctg1[lStartTerm1] == 0L)
                lStartTerm1 += 1L;

        lComputedTerm += 1L;

        lDivider += 2L;

        // Reverse Sign for Arctg series

        bAdd = bAdd ? FALSE : TRUE;

    }

    // Clean Up

    delete [] plArctg1;
    delete [] plArctg2;
    delete [] plArctg3;
}